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1.  SUMMARY 


A  numerical  simulation  capability  has  been  developed  for  the 
analysis  of  fluid  flow  in  liquid-filled  projectiles.  The  unsteady, 
compressible  Navier-Stokes  equations  for  laminar  flow  are  solved  for 
the  primitive  variables  without  recourse  to  linearization  or  simpli¬ 
fication  of  the  equations  of  motion.  The  finite  difference  approxi¬ 
mations  to  the  governing  equations  are  by  choice  either  of  first  or 
second  order  time  accuracy,  second  order  spatial  accuracy  in  the  axial 
and  radial  directions,  and  fourth  order  accuracy  in  the  azimuthal 
direction.  The  method  allows  imposition  of  arbitrary  body  motions 
including  spin  and  precession  and  the  corresponding  boundary  conditions 
are  easily  and  directly  prescribed.  The  finite  difference  equations 
are  solved  using  an  implicit,  approximate  factorization  procedure 
that  permits  the  choice  of  reasonably  large  time  steps  and  avoids 
limitations  based  on  the  magnitude  of  Reynolds  number.  Thus,  the 
numerical  simulation  methodology  considered  provides  a  complete, 
accurate  and  flexible  framework  for  the  computational  analysis  of 
fluid  behavior  in  liquid-filled  projectiles. 

In  the  current  effort,  two  computer  programs  have  been  developed. 
The  cartesian  velocity  components  and  pressure  are  the  dependent 
variables  in  the  first  along  with  a  cartesian  base  coordinate  system. 
This  code  retains  the  desired  features  of  arbitrary  geometry  and  body 
motion.  The  second  computer  program  is  cast  with  cylindrical  velocity 
components  in  cylindrical  coordinates  and  is  currently  (but  not 
inherently)  limited  to  a  cylindrical  geometry.  Both  these  programs 
are  applicable  to  three-dimensional  flow  and  both  codes  can  be  flagged 
to  compute  axi symmetric  flow  efficiently  by  avoiding  calculations 
involving  the  extra  coordinate.  This  axisynmetric  option  has  been 
extensively  exercised  and  the  results  compare  favorably  with  other 
known  numerical  results.  The  fully  three-dimensional  codes  have 
been  exercised  for  fewer  cases  but  enough  to  validate  their  accurate 
applicability. 
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While  it  was  only  required  that  either  a  properly  formulated 
cartesian  variable  or  a  cylindrical  variable  code  be  constructed  that 
can  efficiently  treat  axisymmetric  cases,  both  approaches  have  been 
executed  in  the  present  effort.  The  cartesian  variable  computer 
program  solved  the  governing  equations  in  conservation  form  while  the 
cylindrical  variable  formulation  used  the  non-conservation  form  and 
consequently  somewhat  different  discretizations  were  employed  in  the 
two  codes  as  an  exercise.  The  cylindrical  variable  formulation  gave 
better  results,  but  the  precise  reason  for  the  improvement  is  incon¬ 
clusive  yet. 
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2.  INTRODUCTION 


BACKGROUND 

The  Army  is  currently  developing  liquid  payloads  for  spin- 
stabilized  projectiles.  A  liquid  filler  can  destabilize  the  flight  of 
an  aeroball istically  well  designed  shell.  The  instability  mechanism  is 
a  resonance  between  the  coning  frequency  of  the  projectile  and  a  natural 
vibration  frequency  of  the  spinning  liquid.  Only  simplified  models  have 
been  used  in  the  past  to  study  this  liquid-induced,  projectile  instabi¬ 
lity.  These  include  linear  models  and  axi symmetric  (independent  of 
circumferential  or  azimuthal  position  within  the  fluid)  finite- 

difference  solutions  to  the  governing  equations,  the  Navier-Stokes 
12  3 

equations  The  linear  models  consist  of  closed  form  analytic 

solutions  and/or  numerical  solutions.  The  primary  deficiency  with  the 
available  models  is  that  the  projectile  yaw  must  be  quite  small.  In 
fact,  experiments  to  verify  these  modes  indicate  serious  nonlinear 
effects  at  yaw  amplitudes  well  below  those  common  to  projectiles  . 

Hence,  large  yaw  effects  must  be  incorporated  into  the  modeling  capa¬ 
bilities.  Also,  many  of  the  linear  models  do  not  properly  treat  the 
boundary  conditions  at  the  liquid/solid  interfaces.  Corrections  are 
required  within  these  models,  but  the  corrections  are  only  valid  for 
high  Reynolds  numbers.  The  axi symmetric  Navier-Stokes  code  cannot 
model  the  three-dimensional  disturbances  within  the  liquid  produced 
by  the  yawing  motion  or  the  three-dimensional  response  of  the  liquid. 

The  natural  frequencies  of  the  liquid  are  truly  three-dimensional 
oscillations  and  most  flow  problems  cannot  be  treated  by  an  axi symmet¬ 
ric  code.  A  full  three-dimensional  solution  to  the  Navier-Stokes 
equations  is  required.  Viscous  effects  must  be  maintained  properly 
in  all  coordinate  directions  to  properly  model  the  rotating  fluid. 

The  code  must  also  have  a  high  Reynolds  number  capability.  This 
requires  the  use  of  a  computational  algorithm  which  is  more  complex 
than  a  standard  finite  difference  solution. 
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This  report  describes  the  development  and  testing  of  a  Navier- 

Stokes  solver  that  fits  the  bill.  The  basic  method  was  first  presented 

5 

by  Steger  and  Kutler  as  an  implicit  finite-difference  procedure  to 
solve  the  incompressible  Navier-Stokes  equations  used  to  model  aircraft 
vortex  wake  flows  of  external  aerodynamics.  This  code  was  modified 
extensively  by  Steger^  in  order  to  be  used  for  the  internal  fluid 
mechanics  within  liquid  filled  shells.  These  changes  were  made  to 
enable  the  computer  code  to  accept  arbitrary  body  motion  and  quite 
general  body  geometry,  to  minimize  computer  storage  space,  and  to 
account  for  internal  boundary  conditions.  While  this  code  was 
verified  to  be  stable  for  spin-up  and  spin-down  for  cylindrical 
geometries  and  showed  correct  flow  field  structure  and  development, 
it  was  not  directly  compared  with  other  theories.  Also,  the  computer 
program  required  too  much  computational  time  to  pursue  detailed  spin- 
up  calculations;  having  been  built  for  arbitrary  body  motion  and  shape, 
it  was  not  efficient  for  axi symmetric  spin-up.  It  was  thus  logical 
to  start  from  this  computer  code,  make  the  necessary  modifications  to 
make  it  more  efficient  for  axisymmetric  flows  and  subject  the  method 
to  thorough  testing. 

OBJECTIVES  AND  TASKS 

Towards  the  overall  objective  of  developing  a  numerical  simulation 
capabil  ity  for  liquid-filled  projectiles  that  is  complete,  accurate  and 
flexible,  this  effort  concentrated  on  modifying  an  existing  finite 
difference  method  which  potentially  could  meet  the  requirements  and 
checking  it  thoroughly  against  known  analytic  and  numerical  solutions 
for  simple  cases.  The  method  was  also  to  be  tested  for  limitations 
and  cost  for  a  few  fully  three-dimensional  cases.  These  goals  were 
broken  up  into  the  following  tasks: 

1)  An  existing  computer  code  for  three-dimensional  incompressible  flow 
was  to  be  modified  to  produce  improved  efficiency  and  accuracy. 
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(This  code  was  developed  under  Scientific  Services  Program  Contract 
No.  DAAG29-76-D-0100,  but  it  was  only  superficially  checked  for 
accuracy  and  operational  capabilities)  The  existing  code  was 
formulated  using  cartesian  velocity  components.  These  velocities 
are  not  convenient  for  a  spinning  projectile  where  the  initial  flow 
is  axisymmetric.  In  order  to  improve  the  computational  efficiency, 
the  code  was  to  be  recast  in  cylindrical  coordinates  or  the  axi¬ 
symmetric  variation  was  to  be  subtracted  out. 

2)  The  code  was  also  to  be  set  up  so  that  it  can  easily  be  executed 
in  an  axisymmetric  mode.  This  capability  would  be  necessary  to 
evaluate  questions  of  numerical  accuracy  and  efficiency  and  to 
establish  Reynolds  number  limits. 

3)  The  code  was  to  be  structured  so  that  the  numerical  scheme  was  of 
either  first  or  second  order  time  accuracy.  Solutions  of  second 
order  accuracy  would  allow  increased  efficiency  by  using  larger 
time  steps,  but  they  are  normally  prone  to  more  numerical  instabi¬ 
lities.  The  advantages  and  disadvantages  of  larger  time  steps 
versus  numerical  instabilities  were  to  be  investigated  for 
rotating  fluid  flow  problems. 

4)  Solutions  from  the  axisymmetric  version  of  the  code  were  to  be 
compared  to  other  known  numerical  results. 

5)  Run  times  for  the  three-dimensional  and  axisymmetric  versions  of 
the  code  were  to  be  compared. 

6)  The  axisymmetric  version  was  to  be  exercised  to  define  the 
limitations  of  the  code  for  projectile  spin  rate  and  Reynolds 
number. 

7)  Similar  to  task  6  but  for  fewer  cases,  spin  rate  and  Reynolds 
number  limitations  were  to  be  investigated  for  the  three-dimen¬ 
sional  code  option. 


11 


REPORT  LAYOUT 


The  rest  of  this  report  addresses  these  tasks,  how  they  were 
accomplished  and  the  results  of  the  study.  Section  3  covers  the 
governing  equations  with  cartesian  velocity  components  as  the  dependent 
variables.  The  numerical  method  for  this  formulation  is  described  in 
Section  4.  Results  for  axisymmetric  spinning  flow  obtained  using  the 
above  formulations  is  given  in  Section  5.  The  cylindrical  variable 
formulation  and  results  are  shown  in  Section  6.  The  treatment  in  this 
chapter  is  brief  in  as  much  as  the  primary  thrust  of  this  effort 
was  to  check  out  the  cartesian  variable  formulation.  Three-dimensional 
results  are  covered  in  Section  7  followed  by  concluding  remarks  and 
a  list  of  references  in  Sections  8  and  9. 
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3.  GOVERNING  EQUATIONS  IN  CARTESIAN  VARIABLES 

The  motion  of  a  liquid  fluid  within  a  container  is  governed  by 
the  incompressible  Navier-Stokes  equations. 

Momentum  equations 

3tq  +  3x(e  -  ev)  +  3y(?  -  ?v)  +  3Z(9  -  9V)  *  0  (la) 


Continuity  equation 


3  II  +  3  V  +  3  W  =  0 

X  J 


(lb) 


where 
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The  continuity  equation  has  no  time  derivative.  To  facilitate 
construction  of  an  implicit  ,  approximately  factored  numerical  method, 
the  continuity  equation  is  modified  to  be 


8tp  +  V  +  8yV  +  9zw)  =  9tP* 


TU  •  •  * 

The  quantities  g,  p  and  the  numerical  algorithm  will  be  chosen  to 
let 

9XU  +  3yV  +  3ZW  -4-  0 

(see  next  Section).  Equations  (la)  and  (lc)  may  now  be  combined  to 
yield 

♦  3x(E-Ev)  *  *  az(®-Ev)  =  (8tp*.  0.  0,  0)TranSP°Se 
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In  the  above,  the  subscript  v  was  used  to  denote  viscous  terms. 
But  for  this  and  other  subscripts  that  will  be  clearly  pointed  out,  the 
subscript  notation  will  be  used  interchangeably  with  the  derivative 
operator  notation  9V,  etc.,  to  denote  derivatives.  Thus,  5  f  =  f  . 

X  XX 

In  order  to  permit  arbitrary  body  motion  and  shape,  new  independent 
variables  (but  not  dependent  variables)  are  introduced. 


t  =  t 
S  =  C(x,t) 

n  =  n(x,y,z,t)  (3) 

c  =  c(x,y,z,t) 


These  coordinates  are  defined  in  Fig.  1.  By  restricting  £  to  be 
a  function  of  only  x  and  t,  computer  storage  is  minimized  and  the 
computation  of  metrics  is  simplfied.  Yet,  little  generality  is  lost. 
Using  Eqs.  (3),  the  transformed  equations  become 


9tD  +  9C(E-EV)  +  3n(F-Fv)  +  ac(G-Gy) 


=  O-1 (3tP*,0,0,0)Transpose 


(4a) 


where 


6  =  U/J 

E  =  (etfi  +  cxt  )/J 

F  =  (n^.D  +  nxE  +  nyF  +  nzG)/J 

G  =  (ctD  +  sxE  +  syF  +  szG)/J 

*v  "  ^x^v 

K  =  +  ny^v  +  ^V)/J 

K  =  (^XEV  +  ^yFV  +  4zGy)/J 


(4b) 
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>  * 


SC80-S436 


Fig.  1.  Sketch  of  interior  flow  domain 
with  coordinate  definition 
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and 


E»  = 


CVU-  +  TlvU  + 

x  4  x  n  x  ^ 

VE  +  V„  +  Vc 

Vi  *  Vn  +  V( 


Fv  =  v 


0 

0 

nwu  +  r,u. 

ri  _  u  +  r  u 

y  n  y  c 

G>4 

II 

'z  n  z  ? 

n„v  +  ?wv„ 

V 

n„  v  +  tv 

y  n  y  c 

z  n  z? 

n„w  +  c„w_ 

n,w  +  c  w 

L  y  n  y  cj 

z  n  z  ? 

(5) 


Throughout  the  above  equations,  5y  and  C2  were  assumed  to  be  zero, 
while  the  remaining  metrics  are  given  by 


(cx/j)  =  i/(x5j),  (nx/J)  =  2cy?  -  y^z^> 

•(cx/J)  =  yc2n  "  V 

(ny/j)  =  x^z^ 

(5y/J)  =  -  X^Z, 

(n2/J)  =  -  x^. 

(c2/0)  = 

Ut/J)  =  -  (CX/J)XT 


(nt/J)  =  -(  (nx/J)xx  +  (ny/J)yx  +  (n2/J)zT  ) 
(st/J)  =  -(  (?x/J)xt  +  (?y/J)yT  +  (?2/J)zt  ) 

(  !/J)  =  -  ys2n) 


17 


The  governing  equations  are  kept  in  a  more  manageable  form  if 
contravariant  velocities  are  introduced  (shown  here  in  the  non-scaled 
form) 

u  ■  «t  +  v 

V  =  nt  +  V  +  nyv  +  nzw  (7) 

w  *  't  +  5xu  *  V  +  'zw 

Then,  the  flux  functions  are  simply  written  as 


eu  +  et(p-6) 

ev  +  nt (p-e ) 

ew  +  ct(p-e) 

uU  +  sxp 

,  F  =  J  1 

uV  +  nxP 

a  _  i 

,  G  =  J  1 

uW  +  cxp 

vU 

vV  +  nyP 

vW  +  c  p 

r 

wU 

w V  +  nzP 

wW  +  CZP 

Note  that  even  though  U,  V,  and  W  are  introduced,  the  dependent  variables 
remain  p,  u,  v,  w  -  i.e.,  the  inertial  or  cartesian  velocity  components 
are  retained. 

The  viscous  terms  have  the  form 
Ev  =  (v/J)(V5-7?  1^  95D  +  V£-Vn  1^  3nD  +  v^-v?  1^  3?D) 

(9) 

Fv  =  (v/J)(vn*V£  1^  +  vn-vn  1^  3^  +  vn-vc  1^  3?D) 

Gv  =  (v/J) (vc-vc  1^  3^D  +  vc-vn  1^  3nD  +  vc-v?  1^  3?D) 
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where  the  matrix  1^  is  the  modified  Identity  matrix: 

'0 

in  *  1  1  OO 

1_ 

and  with  the  use  of  orthogonal  coordinates  (when  possible)  simplify  to 
=  (v/J)  ( 

=  (v/J)  (  +  ny  +  nz 

Gv  =  (v/J)  (  C*  +  sj  + 

In  the  above,  the  equations  have  been  cast  in  nondimensional  form 

with 


)  L,  M 

— m  5 

)  I  9  D 
'  — m  n 


)  in  M 

'  -m  c 


(11) 


u/uref  »  V  =  v/uref  ’ 

w  =  w/uref 

x/xref  *  ^  ”  y/xref  ’ 

*  ■  z/2ref 

*  uref^xref 

P  -  P/ufef 

Re  -  v/ ^Xrefuref ^  * 

6  -  6/u2ef 

(12) 


with  the  subscript  ref  denoting  reference  quantities  and  the  ~ 
are  deleted  in  the  governing  equations  for  simplicity.  For  a  spinning 
cylinder,  the  usual  reference  quantities  are  xrg^  =  a,  the  radius  of 
of  the  cylinder  and  u  -  =  aw  where  u  is  the  rotational  speed. 
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BOUNDARY  CONDITIONS 


Boundary  conditions  require  that  the  velocities  U,  V,  W  =  0  at  the 
walls  (no  slip)  if  the  C»n»C  coordinates  are  rigidly  attached  to  the 
shell.  Thus,  along  the  walls. 


u  =  x,  v  =  y  ,  w  =  z  (j 

where  x,  y,  and  z  are  the  velocities  of  a  particle  attached  to  the 
wall.  Alternately,  spin  about  the  axis  can  be  taken  out  of  the  body 
motion  by  setting^ 


V  = 


(14) 


This  latter  approach  has  the  advantage  of  removing  changes  in  the 
velocity  solely  due  to  coordinate  rotation. 


SUBTRACTING  OUT  AXISYMMETRIC  VARIATION 

The  equations  given  above  are  formulated  in  terms  of  cartesian 
components  of  velocity.  Even  for  completely  spun  up  flow  in  a  shell 
that  is  not  coning  or  precessing,  the  cartesian  velocity  components 
possess  circumferential  gradients  as  opposed  to  cylindrical  velocity 
components  which  only  have  zero  gradients.  This  implies  that,  even  to 
compute  the  axi symmetric  flow  after  spin  up,  in  the  present  formulation, 
one  has  to  solve  the  entire  three-dimensional  equations  in  time.  Also, 
even  for  mildly  non-axi symmetric  flows,  the  number  of  grid  points 
needed  to  accurately  represent  the  cartesian  velocity  components  in  the 
circumferential  direction  is  excessive  compared  to  the  number  necessary 
for  cylindrical  components  of  velocity. 

There  are  two  ways  of  alleviating  this  problem.  One  approach  is  to 
use  cylindrical  components  of  velocity  as  the  dependent  variables.  This 
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is  described  in  Section  5.  The  other  approach  is  to  use  the  current 
cartesian  formulation  but  to  subtract  out  the  true  axisymmetric 
variation. 

In  subtracting  the  true  axisymmetric  variation,  the  governing 
equations  would  be  rewritten 


3D  +  3r(E-E  )  +  3  (F-F,  -  (F  -F  ))  +  3  (G-G  )  =  -(H  -H  ) 
t  £  v  n  a  '  v  va''  v'  '  a  va' 


(15) 


where 


0 

0 

0 

0 

H  -H  =  J"1* 

a  av  n 

V[R^(U-£t)cos<j)  -  R<j)n(V-nt)sin«j» 

+  R  (W-cf)cos(})]  -  p  s i n<()/ ( R<()  ) 

+  vJ_1R~2 

V 

-V[R^(U-£^)sin(|)  +  R<t'r|(V-n^.)cos<|) 

+  R?(W-ct)sin«|»]  -  p  cos<f>/(R4>  ) 

w 

In  Eqs.(15)  and  (16),  the  subscript  a  is  used  to  denote  truly 
axisymmetric  terms.  The  letter  R  is  used  to  denote  the  radius  of 

A  A 

the  point  under  consideration.  The  quantity  3  F  is  equal  to  H  . 

~  f|  d  d 

However,  the  term  3  Fa  on  the  left  hand  side  of  the  equation  will  be 

n  q 

differenced  rather  than  exactly  evaluated  and  it  will  cancel  out 
the  numerical  error  of  differencing  3^F  due  to  axisymmetric 
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variation.  Non-axi symmetric  variation  in  n  is  retained  in  the  terms 

^  A  A 

9  F  -  9  F  while  the  correct  axi symmetric  variation  is  held  in  H  .  To 
n  n  a  J  a 

use  this  approach,  metrics  must  be  evaluated  via  the  procedure  explained 
in  Ref.  7  (Eq.  7).  Note  that  the  viscous  terms  also  contribute  to  the 
true  axi symmetric  variation. 

Once  either  of  these  approaches  are  incorporated,  it  becomes  an  easy 
matter  to  build  a  two-dimensional  option  into  the  computer  code  by 
bypassing  the  n  or  azimuthal  inversions.  This  option  becomes  even 
easier  to  achieve  because  the  present  computer  code  comprises  a  single 
index  array  structure.  The  two-dimensional  program  would  allow  the 
calculation  of  high  Reynolds  number  axisymmetric  spin  up  to  proceed 
until  the  three-dimensional  perturbations  of  the  body  motion  need  to 
be  introduced.  It  can  also  be  used  as  an  economic  code  for  optimizing 
numerical  dissipation,  grid  stretching  and  the  use  of  variable  step 
sizes.  The  axisymmetric  code  will  also  be  used  to  compute  solutions 
to  compare  with  other  known  solutions  such  as  those  presented  by 

3 

Kitchens  . 
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4.  NUMERICAL  METHOD  FOR  CARTESIAN  VARIABLE  FORMULATION 


The  numerical  method  described  below  is  presented  in  detail  in  Ref.  5. 
The  approximate  factorization  difference  equations  in  delta  form  (or 
increment  form)  are  given  below  assuming  orthogonal  coordinates  in 
order  to  simplify  writing  the  viscous  terms. 


[_L + 

hJn+1(6?An  - 

ViVe> 

+ 

«iVs 

] 

•n  + 

h Jn+1 ( 6  Bn  - 
n— 

V2-W 

+ 

e.V  A 
i  n  n 

] 

•II  + 

hJn+1(6cCn  - 

V3*«V 

+ 

e.V  A 

1  K  K 

](ffn+1  -  Dn) 

(I  - 

-  pn) 

(17) 

(1  - 

Jn+1/Jn)Dn 

ATJn+1 

( 6^En  +  6  Fn 

+  6^6n  -  1 

-I  Y,I  6  Dn-6 jJ  6  In) 
ri 1  2— m  ri  5  1  3 — m  £  7 

ee[(V? 

V*  +  (VnAn) 

1  +  (Vc> 

2  ] 

Dn 

where 

Yj  *  (v/J)5j 

Y2  =  (S’  +  S’  +  5§  )(v/J) 

y3  =  K  +  ny  +  nz  ^v/J) 

cQ  and  are  0(At) 

h  =  At  for  Euler  time  differencing 

h  =  At/2  for  trapezoidal  time  differencing 
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VrD  =  D.  -  D .  -  ,  A  D  = 

K  J  J-1  n 

V  =  (Dwl  -  Da_1)/(2AC)  . 

VV  5  {*“k+l+“k**Dk+rDk*  "  f“k+t*k-l 


9 


) ( D k~ D k_ i ) } / { 2 ( An ) 2 }  ,  etc. 


and  ,  etc.  have  similar  definitions  save  for  which  is 
fourth  order  accurate,  i.e., 

V  5  <°k+2  -  4°k+l  +  6°k  -  44-l  +  Dk_2)/(12in) 


In  the  above  equation,  all  of  the  metrics  are  kept  at  time  level  n+h 
for  the  case  of  coordinates  moving  with  time.  The  metrics  are 
centrally  differenced  with  those  in  the  azimuthal  direction  being 
fourth  order  accurate  and  the  others  being  second  order  accurate. 


Also  in  the  above  equation  describing  the  implicit  approximate 
factorization  procedure,  the  equations  have  been  locally  linearized 


in  time  and  the 

Jacobian  matrices  for 

the  inviscid  terms  are  obtained 

from 

L0 

4 

4 

L3 

AAA 

L1 

Q+LjU 

L2U 

L3u 

A,  B,  C 

- 

L2 

4V 

Q+L2v 

L3V 

(18) 

.4 

4W 

L2w 

Q+L3w  ! 

where  Q  =  Lq 

+  L.u  + 

^2v  +  ^3W 

and  Lq  = 

4 

=  4 

»  ^2  5y» 

L3  =  5 

A 

to  obtain  Matrix  A,  etc. 
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Equivalent  matrices  for  the  viscous  terms  may  be  derived  easily  but  are 
avoided  by  using  orthogonal  coordinates  when  possible.  The  terms 
multiplied  by  the  coefficients  eg  and  are  the  explicitly  and 
implicitly  added  smoothing  terms  which  help  stabilize  the  calculations. 

SATISFYING  THE  CONTINUITY  EQUATION 

It  can  be  shown  that  when  @  =  0(1/(at)),  the  equation  (1)  along 
with  the  definitions  (2)  will  satisfy  the  continuity  equation  upto  first 
order  accuracy.  For  better  satisfaction  of  the  continuity  equation,  the 
finite  difference  equations  (14)  are  solved  atleast  twice  every  time 
step.  In  the  first  iteration,  p  is  set  to  be  =  pn  and  6  is  chosen 
to  be  of  the  order  of  1/(at).  For  subsequent  iterations,  p*  is  set 
to  the  value  of  pn  ^  predicted  in  the  previous  iteration  and  6  is 
set  to  unity.  This  strategy  serves  the  double  purpose  of  satisfying 
the  continuity  equation  better  (when  p*  -  pn  =  pn+1  -  pn,  the 
continuity  equation  is  satisfied  exactly)  and  removing  errors  due  to 
large  values  of  $  contaminating  the  momentum  equations  due  to 
approximate  factorization.  For  more  details,  see  Appendix  A  of 
Ref.  5. 
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5.  RESULTS  FOR  AXISYMMETRIC  SPIN-UP 
USING  CARTESIAN  VARIABLE  FORMULATION 

Several  axisymmetric  test  cases  involving  pure  spin  up  from  rest 

were  analyzed  in  order  to  verify  the  method  described  in  the  previous 

Sections.  Results  were  compared  with  solutions  obtained  using  the 

3 

code  developed  by  Kitchens  .  These  latter  solutions  were  provided  by 
Sedney  and  Gerber  of  BRL  (private  communication).  This  Section  is 
devoted  to  a  description  of  the  results  obtained  using  the  cartesian 
variable  formulation.  Results  using  the  cylindrical  variable  formu¬ 
lation  will  be  presented  in  the  next  Section  after  presenting  the 
governing  equations. 

GRID  CLUSTERING 

Near  the  end  and  side  walls  of  the  cylindrical  container,  adequate 
flow  field  resolution  is  especially  important  to  accurately  account 
for  the  predominantly  viscous  effects  in  the  boundary  layers.  To  achieve 
this  resolution,  the  same  formulae  used  by  Kitchens  to  cluster  more 
grid  points  near  the  end  walls  and  side  wall  were  employed.  For  example, 
in  the  radial  direction,  the  radial  position  of  a  grid  point  corres¬ 
ponding  to  a  uniformly  spaced  coordinate  c  'was  defined  by 

{(b+l)/(b-l)}C  -  1 

r  =  b  -  (19) 

{(b+l)/(b-l)}C  +  1 

(the  parameter  b  is  the  control  to  achieve  desired  clustering) 


where  c  and  r  range  from  zero  to  one.  The  axial  direction  was 
similarly  treated. 
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EVOLUTION  OF  SPIN-UP  FLOW 

Results  are  now  presented  to  illustrate  the  evolution  of  the  flow 
field  in  a  cylindrical  spinning  cavity  spun  up  impulsively  from  rest. 

The  first  case  considered  was  a  cylinder  of  aspect  ratio  unity  (aspect 
ratio  is  the  ration  of  cylinder  half  height  to  its  radius)  with  a  fluid 
of  Reynolds  number  10,000  based  on  the  cylinder  radius  and  rotational 
speed  (Reynolds  number  =  u>a2/v;  see  Eq.  (12)  for  notation).  The  finite 
difference  grid  employed  for  this  case  is  shown  in  Figure  2  for  half 
the  cylinder  (31  x  31  points).  In  this  and  the  following  figures,  the 
axis  labelled  R  denotes  the  radial  coordinate  and  the  axis  labelled 
Z  corresponds  to  the  axial  coordinate  (contrary  to  the  convention  used 
in  Fig.  1). 

The  in-plane  (radial  and  axial)  velocity  field  is  plotted  in  Fig.  3a 
from  the  solution  after  125  steps  at  a  non-dimensional  time  step  of  0.25 
(non-dimensional  time  t  =  31.25).  The  maximum  in-plane  velocity 
component  was  'v  0.15  and  a  length  segment  of  0.075  units  was  chosen  to 
represent  this  magnitude  of  velocity.  All  velocity  vectors  were  scaled 
with  respect  to  this  value.  Arrows  were  drawn  centered  at  each  grid 
point  with  their  length  representing  the  magnitude  of  the  velocity  there 
and  the  orientation  depicting  direction  of  velocity.  While  there  is 
very  little  action  in  most  of  the  cylinder,  the  fluid  particles  close 
to  the  end  walls  (axial  ends)  are  clearly  seen  to  be  accelerated 
outward.  Two  close  ups  of  the  same  plot  focusing  on  one  corner  are 
shown  in  Figs.  3b  and  3c  for  better  clarity.  The  same  length  scales 
are  used  in  all  velocity  vector  plots  (a  magnitude  of  0.15  is  scaled 
to  a  length  of  0.075)  for  this  case  to  enable  easy,  quantitative 
visualization  of  the  flow  processes.  Figures  3d-3h  show  the  decay  in 
the  magnitude  of  the  in-plane  velocity  field  as  the  fluid  spins  up  and 
approaches  the  steady  state  of  solid  body  rotation. 

It  is  also  interesting  to  trace  the  evolution  of  the  azimuthal 
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COMPUTATIONAL  GRID 
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Fig.  2.  Computational  Grid  for  Reynolds  number  of 

10,000. 
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IN-PLANE  VELOCITY  FIELD 


IN-PLANE  VELOCITY  FIELD 


Fig.  3a.  Time  =  31.25 


Fig.  3b.  Close  up  of  Fig.  3a 


IN-PLANE  VELOCITY  FIELD 


IN-PLANE  VELOCITY  FIELD 


Fig.  3c.  Close  up  of  Fig.  3b 


Fig.  3d.  Time  =  62.5 


Fig.  3. 


In-plane  velocity  field 
Reynolds  Number  =  10,000. 


(page  1  of  2) 
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IN-PLANE  UELOCITY  FIELD 


IN-PLAHE  UELOCITY  FIELD 


Fig.  3e.  Time  =  93.75 


Fig.  3f ,  Time  =  156.25 


IN-PLANE  UELOCITY  EIELD 


R 

Fig.  3g.  Time  =  218.75 


IN-PLANE  UELOCITY  FIELD 


Fig.  3h.  Time  =  281.25 


Fig.  3.  (concluded) 


30 


velocity  distribution.  The  sequence  of  Figs.  4a  -  4f  follow  the 
progression  of  this  velocity  component  as  the  flow  field  develops 
towards  solid  body  rotation.  The  minimum  contour  level  shown  is  0.05 
and  the  maximum  value  is  1.0  with  increments  between  levels  of  0.05. 

A  stream-function  was  defined  by 

9(4<)/9(radial  coordinate)  =  (axial  velocity) 

(20) 

8(if<)/9(axial  coordinate)  =  -(radial  velocity) 


and  thus 


V2\p  =  vorticity  (21) 

Streamline  plots  were  drawn  at  a  sequence  of  time  steps  and  are 
displayed  in  Figs.  5a  -  5c.  Tangents  drawn  to  these  contours  must 
be  parallel  to  the  in-plane  velocity  at  that  point  and  it  is 
interesting  to  observe  the  vortex  core  move  from  the  corner  between 
end  and  side  walls  radially  inward  as  the  flow-field  develops. 

Similar  sets  of  plots  are  shown  in  Figures  6-9  for  a  Reynolds 
number  of  1,000.  Figure  6  displays  the  computational  grid  of  points 
(31  x  31).  A  less  fine  clustering  of  points  is  evident  compared  to 
the  grid  for  Reynolds  number  10,000.  The  development  of  the  in-plane 
velocity  field  and  the  azimuthal  velocity  contours  are  portrayed  in 
Figs.  7a  -  7c  and  Figs.  8a  -  8b  respectively  after  250  and  500  steps 
at  a  step  size  of  0.25.  Streamlines  are  drawn  after  250  steps  in 
Fig.  9.  This  lower  Reynolds  number  case  is  almost  fully  spun-up 
to  the  solid  body  state  of  rotation  after  only  500  steps.  This  is 
in  contrast  to  the  first  case  which  converges  much  more  slowly  to 
the  asymptotic  steady  state. 
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AZIMUTHAL  VELOCITY  CONTOURS 


AZIMUTHAL  VELOCITY  CONTOURS 
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Fig.  4a.  Time  =  31.25 


Fig.  4b.  Time  =  62.5 
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Fig,  4c.  Time  =  93,75 


Fig.  4d,  Time  =  156.25 


Fig.  4.  Azimuthal  velocity  plots  for 
Reynolds  number  10,000 
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AZIMUTHAL  VELOCITY  COMTOURS 


AZIMUTHAL  VELOCITY  COMTOURS 


Fig.  4e.  Time  =  218,75 


Fig.  4f .  Time  =  281.25 


Fig.  4.  (concluded) 


STREAMLINE  CONTOURS 
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Fig.  5a.  Time  =  31.25 


Fig.  5b,  Time  =  62.5 


Fig.  5.  Streamline  contour  plots 

for  Reynolds  number  10,000  (page  1  of  2) 
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STHEAALIFC  CONTOURS 
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Fig.  5c,  Time  =  93.75 


Fig.  5.  (concluded) 
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COMPUTATIONAL  GRID 
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Reynolds  number  1,000 
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IN-PLANE  VELOCITY  FIELD 
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Fig.  7a.  Time  =  62.5 


Fig.  7bt  Close  up  of  Fig. 
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AZIHUTHAl  UEIOCITY  CONTOURS 


AZINUTHAl  UEIOCITV  CONTOURS 


Fig.  8a.  Time  =62.5 


Fig.  8b.  Time  =  125.0 


Fig.  8.  Azimuthal  velocity  plots 
for  Reynolds  number  1,000 
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Fig.  9.  Streamline  contours  at  Time  =  62.5 
for  Reynolds  number  1,000. 
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COMPARISON  WITH  RESULTS  FROM  KITCHENS'  CODE 


A  series  of  comparisons  are  now  made  between  the  results  from  the 

O 

method  under  evaluation  and  the  code  developed  by  Kitchens  .  For  Reynolds 
number  10,000  Figs.  10a  -  lOd  compare  the  two  solutions  at  various  axial 
locations  by  plotting  the  radial  variation  of  azimuthal  velocity.  For 
these  plots,  t  =  62.5  with  At  =  0.10.  The  aspect  ratio  of  the  cylinder 
was  unity.  First  order  time  accuracy  was  used  along  with  two  iterations 
every  time  step  to  correct  the  error  in  satisfying  the  continuity  equation. 
The  two  solutions  match  well,  even  if  not  exactly. 

A  time  history  of  the  azimuthal  velocity  at  one  point  (radial 
location  =  0.9,  axial  location  =  1.0)  is  shown  in  Fig.  11.  At  time 
t  =  62.5,  points  from  this  figure  correspond  to  points  from  Fig.  lOd. 
Looking  at  Fig.  11,  independently  of  Figs.  10,  it  appears  like  the 
present  solution  is  lagging  behind  Kitchens'  solution.  But  ,  comparing 
Figs.  10  and  11,  it  is  seen  that  the  time  history  trends  match  well 
if  not  superbly. 

Similar  comparisons  are  made  for  Reynolds  number  1,000  in  Figs. 

12  and  13  at  two  non-dimensional  times  of  62.5  and  125  with  a  non- 
dimensional  time  step  of  0.25.  Greater  smoothing  was  required  for 
this  larger  time  step  and  compromised  the  solution  accuracy  at  t  =  62.5. 

At  the  later  time  of  125.0,  the  fluid  has  almost  fully  spun  up  and  the 
discrepancy  between  the  two  solutions  is  very  minor. 
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Fig.  10.  Comparison  of  azimuthal  velocity  profiles  for 
Reynolds  number  10,000  and  Time  =  62.5 
(page  1  of  2) 
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AZIMUTHAL  VELOCITY 


SC81-15658 


NON  DIMENSIONAL  TIME 


Fig.  11.  Time  history  of  azimuthal  velocity 
at  axial  location  of  1.0  and 
radial  location  of  0.9. 
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Fig.  13.  Comparison  of  azimuthal  velocity  profiles 
for  Reynolds  number  1,000  at  Time  =  125.0 


6.  CYLINDRICAL  VARIABLE  FORMULATION 


The  governing  Navier-Stokes  momentum  and  continuity  equations 
may  be  written  in  non-divergence  form  in  cylindrical  variables  as 

Continuity  equation 

3ru  +  (l/r)30v  +  3zw  +  u/r  =  0  (22a) 

Momentum  equations 

3tu+u3ru+(v/r)3eu+w3zu-v2/r+  3fp  =  v(V2u-(2/r2)3Qv-u/r2)  (22b) 

3tv+u3rv+(v/r)3gV+w3zv+uv/r+(l/r)3Qp  =  v(V2v+(2/r2)30u-v/r2) 

3tw+u3rw+(v/r)3gW+w3zw  +  3ZP  =  v(V2w  ) 

where,  .  . .  „  , 

V2  =  32  +  (l/r)3r  +  (l/r2)32  +  3| 

r  =  radial  u  =  radial  velocity 

6  =  azimuthal  coordinate  ,  v  =  azimuthal  velocity 

z  =  axial  coordinate  ,  w  =  axial  velocity 


( 


Once  again,  the  continuity  equation  is  modified  by  the  addition 
of  transient  terms  for  pressure 


3tP  +  B(3fu  +  (l/r)30v  +  3zw  +  u/r)  =  3tp 


(22c) 
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These  equations  were  coded  into  an  implicit  factorization  procedure 
modified  from  the  method  presented  in  Section  4.  The  non-conservation 
form  allowed  easy  use  of  upwind  differencing  of  the  convective  terms, 
thus  eliminating  the  need  for  artificial  smoothing  of  the  type  used  in 
Section  4  and  avoiding  dissipation  in  the  continuity  equation  altogether. 
The  results  that  will  now  be  presented  for  axi symmetric  spin  up  show 
a  marked  improvement  over  the  results  using  the  cartesian  variables. 
However,  it  is  too  soon  to  identify  specific  differences  between  the 
two  methods  that  cause  this  improvement. 

RESULTS 

The  parameters  of  the  first  case  discussed  are:  Reynolds  number 
=  1,000,  aspect  ratio  *  1.0.  Comparisons  between  solutions  from  Kitchens' 
method  and  the  current  cylindrical  variable  formulation  are  made  in 
Figs.  14  -  15  at  two  axial  locations  (  z  =  0.0662  ,  z  =  1.0)  at  two 
time  levels  (250  steps  and  500  steps  with  step  size  =  0.25).  A  good 
match  is  observed,  better  than  the  agreement  for  the  results  using  the 
cartesian  variable  formulation. 

The  Reynolds  number  is  10,000  and  the  aspect  ratio  is  unity  for  the 
next  case  considered.  Results  are  compared  at  two  axial  locations 
(z=0.0324  and  1.0)  and  three  time  levels  (250,  500,  1000  steps  with 
step  size  =  0.25)  in  Figs.  16,  17,  and  18.  Once  again,  better 
agreement  between  the  two  solutions  is  seen  compared  to  the  results 
of  the  previous  section.  The  improvement  is  in  spite  of  using  a 
larger  time  step  than  used  for  the  cartesian  variable  formulation. 

The  effect  of  different  values  of  the  parameter  6  and  the 
time  accuracy  are  plotted  in  Figs.  19  and  20.  Both  6  and  the 
time  accuracy  are  seen  to  have  negligible  influence  on  the  results. 

Finally,  as  a  check  to  determine  any  Reynolds  number  limits 
(spin  rate  limits  are  related  to  Reynolds  number  limits)  on  the 
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Fig.  15.  Comparison  of  azimuthal  velocity  profiles 
for  Reynolds  number  1,000  at  Time  =  125.0 
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for  Reynolds  number  10,000  at  Time  =  62.5 
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Fig.  17.  Comparison  of  azimuthal  velocity  profiles 
for  Reynolds  number  10,000  at  Time  =  125.0 
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Fig.  19.  Effect  of  the  parameter  6  shown  for 
Reynolds  number  10,000  at  Time  =62.5 
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Fig.  20.  Effect  of  time  accuracy. 
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computational  method,  results  were  obtained  for  500  steps  (step  size  = 
0.25)  for  a  Reynolds  number  of  a  million.  The  computational  grid  for 
this  case  (41  x  41  points)  is  shown  in  Fig.  21a  and  azimuthal  velocity 
profiles  are  shown  after  250  and  500  steps  in  Fig.  21b.  No  instabi¬ 
lities  were  observed,  lending  credence  to  the  prediction  that  the 
implicit  scheme  under  consideration  would  not  have  any  inherent 
limitations  based  on  the  magnitude  of  Reynolds  number. 
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7,  THREE-DIMENSIONAL  RESULTS 


Both  the  cartesian  and  cylindrical  variable  computer  programs  have 
been  tested  in  the  three-dimensional  mode  of  operation  (circumferential 

i 

direction  is  included  in  the  calculations),  but  for  computing  only 
axisymmetric  flows.  The  axisymmetric  results  presented  in  the  earlier 
sections  could  all  be  reproduced  but  at  greater  computational  cost.  No 
instabilities  unique  to  the  three-dimensional  calculations  were  observed 
and  the  methods  were  stable  for  all  Reynolds  numbers  considered  (100; 

1000;  10,000;  1,000,000). 

A  discussion  of  the  computational  speed  would  conclude  this  section. 
The  nucleus  or  primary  building  block  which  must  be  known  to  compute 
computational  time  for  any  particular  case  is  the  time  for  one  grid 
point  for  one  iteration  to  be  denoted  by  a.  For  the  computer  programs 
under  consideration,  a  =  0,00018.  Knowing  a  and  the  number  of  times 
this  unit  is  repeated  in  a  calculation,  the  total  time  required  may  be 
computed.  For  implicit,  approximately  factored  schemes,  the  calcula¬ 
tions  are  broken  up  into  the  solution  of  one-dimensional  block  tridiagonal 
(4X4  blocks)  system  of  linear  algebraic  equations.  Considering  for 
example,  each  line  of  points  along  which  only  the  index  J  varies,  the 
time  for  one  inversion  would  be  proportional  to  JMAX,  the  total  number 
of  J  points.  There  would  be  KMAX  *  LMAX  such  lines  to  be  treated  in 
the  first  factor  of  the  approximate  factorization.  Thus,  for  calculating 
each  factored  segment  (see  Eq.  17)  in  the  implicit  method,  the  computa¬ 
tional  time  would  be  proportional  to  the  total  number  of  points 
(=  JMAX*KMAX*LMAX) .  The  next  point  to  remember  is  that  periodic  block 
tridiagonal  inversions  take  about  1.7  times  longer  than  non-periodic 
inversions.  Thus  the  circumferential  direction  will  be  proportionately 
more  expensive.  Now,  neglecting  the  time  involved  in  computing  the 
right  hand  side  of  Eq.  (17)  in  comparison  with  the  block  tridiagonal 
inversion  process,  the  formula  that  can  be  used  to  calculate  an  approximate 
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computational  time  for  any  case  is  given  by 


T  =  (JMAX  *  KMAX  *  LMAX) 

*(2  +  I K* 1.7)  (23) 

*  NUMBER  OF  3  ITERATIONS 

*  a 

where 

T  =  total  computational  time  (CPU  time) 

a  =  CPU  time  for  1  grid  point  for  1  iteration  =  0,00018 

(for  CDC  7600) 

1^  =  0  for  two-dimensional  calculations 

(no  circumferential  direction  included  in  calculations) 

=  1  for  three-dimensional  calculations 

(circumferential  direction  included  In  calculations 
with  periodic  boundary  conditions) 

Thus  for  example,  a  two-dimensional  calculation  for  250  time  steps 
with  a  grid  of  31  X  31  points  and  two  .  B  Iterations  per  step  would 
take  173  CPU  seconds  on  a  CDC  7600.  A  21  X  21  X  20(c1rcumferential ) 
grid  would  need  1175  CPU  seconds  for  100  time  steps  with  two  B  iterations 
each  step.  In  contrast,  the  code  developed  by  Kitchens  would  only 
take  37  seconds  for  the  first  example  (31  X  31  points,  250  time  steps). 
However,  the  Kitchens  technique,  being  based  on  an  axi symmetric  stream- 
function,  cannot  be  directly  extended  to  three-dimensional  problems. 


56 


8.  CONCLUDING  REMARKS 


This  report  describes  the  successful  development  of  a  numerical 
simulation  capability  for  liquid  filled  projectiles  applicable  to  three- 
dimensional  flow.  The  emphasis,  in  this  effort,  has  been  on  proving  the 
method  by  comparison  with  other  known  results  for  axi symmetric  flow  and 
testing  the  scheme  for  limitations  and  cost  when  operated  in  the  fully 
three-dimensional  mode.  In  terms  of  the  task  breakdown  given  in  Section 
2  (Introduction),  Tasks  1-3  have  been  covered  in  Sections  3  and  4,  Tasks 
4  and  6  have  been  discussed  in  Sections  5  and  6,  and  Tasks  5  and  7  have 
been  dealt  with  in  Section  7.  Two  formulations  have  been  developed,  one 
based  on  the  cartesian  variables  and  the  other  based  on  cylindrical 
variables.  The  cylindrical  variable  formulation  compares  better  with 
other  numerical  results  for  axi symmetric  spin-up.  Both  formulations 
have  proved  themselves  to  be  viable  methods  to  be  exploited  in  future 
work  to  compute  complex  three-dimensional  flows  with  both  spin  and 
precession. 
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Madison,  IN  47251 

1  Commander 

US  Army  Research  Office 
ATTN:  Dr.  R.E.  Singleton 
P.0.  Box  12211 
Research  Triangle  Park, 

NC  27709 

1  AGARD-NATO 

ATTN:  R.H.  Korkegi 
APO  New  York  09777 

1  Director 

US  Army  TRADOC  Systems 
Analysis  Activity 
ATTN:  ATAA-SL,  Tech  Lib 
White  Sands  Missile  Range 
NM  88002 

3  Commander 

Naval  Air  Systems  Command 
ATTN:  AIR-604 
Washington,  DC  20360 

2  Commander 

David  W.  Taylor  Naval  Ship 
Research  &  Development  Center 
ATTN:  H. J .  Lugt,  Code  1802 

S.  de  los  Santos,  Head, 
High  Speed  Aero  Division 
Bethesda,  MD  20084 


1  Commander 

Naval  Surface  Weapons  Center 
ATTN:  DX-21 ,  Lib  Br 
Dahlgren,  VA  22448 

4  Commander 

Naval  Surface  Weapons  Center 
Applied  Aerodynamics  Division 
ATTN:  K.R.  Enkenhus 
M.  Ciment 

A. E.  Winklemann 
W.C.  Ragsdale 

Silver  Spring,  MD  20910 

1  AFATL  (DLDL,  Dr.  D.C. Daniel) 
Eglin  AFB,  FL  32542 

2  AFFDL  (W.L.  Hankey;  J.S.  Shang) 
Wright-Patterson  AFB,  OH  45433 

5  Director 

National  Aeronautics  and  Space 
Admi nistration 
ATTN:  D.R.  Chapman 
J.  Rakich 
W.C.  Rose 

B.  Wick 
P.  Kutler 

Ames  Research  Center 
Moffett  Field,  CA  94035 

4  Director 

National  Aeronautics  and  Space 
Administration 
ATTN:  E.  Price 
J.  South 
J.R.  Sterrett 
Tech  Library 
Langley  Research  Center 
Langley  Station 
Hampton,  VA  23365 

1  Director 

National  Aeronautics  and  Space 
Administration 

STAD-Theoretical  Aero.  Branch 
ATTN:  Dr.  M.  D.  Salas 
Hampton,  VA  23365 


60 


No.  of 
Copies 


DISTRIBUTION  LIST 


No.  of 

Organi zati on  Copies  Organi zati on 


1  Di rector 

National  Aeronautics  and  Space 
Administration 
Lewis  Research  Center 
ATTN:  MS  60-3,  Tech  Lib 
21000  Brookpark  Road 
Cleveland,  OH  44135 

2  Di rector 

National  Aeronautics  and 
Space  Administration 
Marshall  Space  Flight  Center 
ATTN:  A.R.  Felix,  Chief 
S&E-AERO-AE 
Dr.  W.W.  Fowl  is 
Huntsville,  AL  35812 

2  Director 

Jet  Propulsion  Laboratory 
ATTN:  L.M.  Mach 

Tech  Library 
4800  Oak  Grove  Drive 
Pasadena,  CA  91103 

3  Arnold  Research  Org.,  Inc. 
ATTN:  J.D.  Whitfield 

R.K.  Matthews 
J.C.  Adams 

Arnold  AFB,  TN  37389 

3  Aerospace  Corporation 
ATTN:  H.  Mirels 

R.L.  Varwig 
Aerophysics  Lab. 

P.0.  Box  92957 

Los  Angeles,  CA  90009 

1  AVCO  Systems  Division 
ATTN:  B.  Reeves 
201  Lowell  Street 
Wilmington,  MA  01887 

3  Boeing  Commercial  Airplane 
Company 

ATTN:  R.  A.  Day,  Org  B-8120 

P.E.  Rubbert,  MS  3N-19 
J.D.  McLean,  MS-3N-19 
Seattle,  WA  98124 


3  Cal  span  Corporation 
ATTN:  A.  Ritter 
G.  Homicz 
W.  Rae 
P.0.  Box  400 
Buffalo,  NY  14225 

1  General  Dynamics 

ATTN:  Research  Lib  2246 

P.0.  Box  748 

Fort  Worth,  TX  76101 

1  General  Electric  Company,  RESD 
ATTN:  W.  J.  East 

3198  Chestnut  Street 
Philadelphia,  PA  19101 

2  Grumman  Aerospace  Corporation 
ATTN:  R.E.  Melnik 

L.G.  Kaufman 
Bethpage,  NY  11714 

2  Lockheed-Georgia  Company 
ATTN:  B.H.  Little,  Jr. 

G. A.  Pounds 

Dept  72074,  Zone  403 
86  South  Cobb  Drive 
Marietta,  GA  30062 

1  Lockheed  Missiles  and  Space 

Company 

ATTN:  Tech  Info  Center 
3251  Hanover  Street 
Palo  Alto,  CA  94304 

3  Martin-Marietta  Laboratories 
ATTN:  S.H.  Maslen 

S.C.  Traugott 

H.  Obremski 
1450  S.  Rolling  Road 
Baltimore,  MD  21227 

2  McDonnell  Douglas  Astronautics 

Corporation 
ATTN:  J.  Xerikos 
H.  Tang 

5301  Bolsa  Avenue 
Huntington  Beach,  CA  92647 
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2  McDonnel 1 -Dougl as  Corporation 
Douglas  Aircraft  Company 
ATTN:  T.  Cebeci 

K.  Stewartson 
3855  Lakewood  Boulevard 
Long  Beach,  CA  90801 

2  Rockwell  International  Science 
Center 

ATTN:  Dr.  V.  Shankar 
Dr.  N.  Malmuth 
1049  Camino  Dos  Rios 
Thousand  Oaks,  CA  91360 

2  Sandia  Laboratories 
ATTN:  F.G.  Blottner 
Tech  Lib. 

Albuquerque,  NM  87115 

1  Union  Carbide  Nuclear  Division 
ATTN:  Dr.  J.E.  Park,  MS  53 

Bldg.  K-1007 
P.0.  Box  P 

Oak  Ridge,  TN  37830 

2  United  Aircraft  Corporation 
Research  Laboratories 
ATTN:  M.J.  Werle 

Library 

East  Hartford,  CT  06108 

1  LVT  Aerospace  Corp. 

Vought  Systems  Division 
ATTN:  J.M.  Cooksey,  Chief, 

Gas  Dynamics  Lab, 
2-53700 
P.0.  Box  5907 
Dallas,  TX  75222 

1  Arizona  State  University 
Department  of  Mechanical  and 
Energy  Systems  Engineering 
ATTN:  G.P.  Neitzel 
Tempe,  AZ  85281 


3  California  Institute  of 
Technol ogy 
ATTN:  Tech  Library 
H.B.  Keller 

Mathematics  Dept. 

D.  Coles 

Aeronautics  Dept. 
Pasadena,  CA  91102 

1  Cornell  University 

Graduate  School  of  Aero  Engr 
ATTN:  Library 
Ithaca,  NY  14850 


1  Illinois  Institute  of  Tech 
H.  M.  Nagib 
3300  South  Federal 
Chicago,  IL  60616 

1  The  Johns  Hopkins  University 
Department  of  Mechanics  and 
Materials  Science 
ATTN:  S.  Corrsin 
Baltimore,  MD  21218 

4  The  Johns  Hopkins  University 
Applied  Physics  Laboratory 
ATTN:  Dr.  R.D.  Whiting 
Dr.  D.A.  Hurdif 
Dr.  R.S.  Hirsh 
Mr.  E.R.  Bohn 
Johns  Hopkins  Road 
Laurel,  MD  20707 

1  University  of  California 
Lawrence  Livermore  National 
Laboratory 
ATTN:  Dr.  Viecelli 
Livermore,  CA  94550 

3  Massachusetts  Institute  of 
Technology 
ATTN:  E.  Covert 

H.  Greenspan 
Tech  Lib 

77  Massachusetts  Avenue 
Cambridge,  MA  02139 
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2  North  Carolina  State  Univ 
Mechanical  and  Aerospace 
Engineering  Department 
ATTN:  F.F.  DeJarnette 
J.C.  Williams 
Raleigh,  NC  27607 

1  Notre  Dame  University 
Department  of  Aero  Engr 
ATTN:  T.J.  Mueller 
South  Bend,  IN  46556 

2  Ohio  State  University 
Dept  of  Aeronautical  and 

Astronautical  Engineering 
ATTN:  S.L.  Petrie 

O.R.  Burggraf 
Columbus,  OH  43210 

2  Polytechic  Institute  of 

New  York 

ATTN:  G.  Moretti 
S.G.  Rubin 
Route  110 

Farmingdale,  NY  11735 

3  Princeton  University 

James  Forrestal  Research  Ctr 
Gas  Dynamics  Laboratory 
ATTN:  S.M.  Bogdonoff 
S.I.  Cheng 
Tech  Library 
Princeton,  NJ  08540 

1  Purdue  University 

Thermal  Science  &  Prop  Ctr 
ATTN:  Tech  Library 
W.  Lafayette,  IN  47907 

1  Rensselaer  Polytechnic 
Institute 

Department  of  Math  Sciences 
ATTN:  R.C.  Diprima 
Troy,  NY  12181 


1  San  Diego  State  University 
Department  of  Aerospace  Engr 
and  Engineering  Mechanics 
College  of  Engineering 
ATTN:  K.C.  Wang 
San  Diego,  CA  92182 

1  Southern  Methodist  University 
Department  of  Civil  and 
Mechanical  Engineering 
ATTN:  R.L.  Simpson 
Dallas,  TX  75275 

1  Southwest  Research  Institute 
Applied  Mechanics  Reviews 
8500  Culebra  Road 
San  Antonio,  TX  78228 

6  Stanford  University 

Dept  of  Aeronautics/ Astronautics 
ATTN:  Dr.  J.L.  Steger  (3  cys) 
Dr.  S.  Chakravarthy 
(3  cys) 

Stanford,  CA  94305 

1  Texas  A&M  University 
College  of  Engineering 
ATTN:  R.H.  Page 
College  Station,  TX  77843 

1  University  of  California  - 
Berkel ey 

Department  of  Aerospace 
Engineering 
ATTN:  M.  Holt 
Berkeley,  CA  94720 

1  University  of  California  - 
Davis 

ATTN:  H.A.  Dwyer 
Davis,  CA  95616 
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2  University  of  California  - 
San  Diego 

Department  of  Aerospace 
Engineering  and  Mechanical 
Engineering  Sciences 
ATTN:  P.  Libby 

Tech  Library 
La  Jolla,  CA  92037 

1  University  of  Cincinnati 
Department  of  Aerospace 
Engineering 
ATTN:  R.T.  Davis 
Cincinnati,  OH  45221 

1  University  of  Colorado 

Department  of  Astro-Geophysics 
ATTN:  E.R.  Benton 
Boulder,  CO  80302 

1  University  of  Hawaii 
Dept  of  Ocean  Engineering 
ATTN:  G.  Venezian 
Honolulu,  HI  96822 

2  University  of  Maryland 
ATTN:  W.  Melnik 

J.D.  Anderson 
College  Park,  MD  20742 

2  University  of  Michigan 
Department  of  Aeronautical 
Engineering 
ATTN:  W.W.  Wilmarth 
Tech  Library 

East  Engineering  Building 
Ann  Arbor,  MI  48104 

1  University  of  Santa  Clara 
Department  of  Physics 
ATTN:  R.  Greeley 
Santa  Clara,  CA  95053 
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3  University  of  Southern 
Cal i fornia 

Department  of  Aerospace 
Engi neering 
ATTN:  T.  Maxworthy 
P.  Weidman 
L.G.  Redekopp 
Los  Angeles,  CA  90007 

1  University  of  Texas 
Department  of  Aerospace 
Engi neeri ng 

ATTN:  J.C.  Westkaemper 
Austin,  TX  78712 

1  University  of  Virginia 
Department  of  Aerospace 
Engineering  &  Engineering 
Physics 

ATTN:  I.D.  Jacobson 
Charlottesville,  VA  22904 

1  University  of  Virginia  Research 
Laboratories  for  the 
Engineering  Sciences 
ATTN:  Prof.  H.  G.  Wood 
P.0.  Box  3366 
University  Station 
Charlottesville,  VA  22903 

1  University  of  Washington 
Department  of  Mechanical 
Engi neeri ng 
ATTN:  Tech  Library 
Seattle,  WA  98105 

1  University  of  Wyoming 
ATTN:  D.L.  Boyer 
University  Station 
Laramie,  WY  82071 

1  Virginia  Polytechnic  Institute 
and  State  University 
Department  of  Aerospace 
Engineering 
ATTN:  Tech  Library 
Blacksburg,  VA  24061 


64 


No.  of 
Copies 


DISTRIBUTION  LIST 


Organization 

1  Woods  Hole  Oceanographic 
Institute 

ATTN:  J.A.  Whitehead 
Woods  Hole,  MA  02543 


Aberdeen  Proving  Ground 

Director,  USAMSAA 
ATTN:  DRXSY-D 

DRXSY-MP,  H.  Cohen 

Commander,  USATECOM 
ATTN:  DRSTE-TO-F 

Director,  USACSL,  Bldg.  E3330,  EA 
ATTN:  Munitions  Division 
W.C.  Dee 

Director,  USACSL,  Bldg.  E3516,  EA 
ATTN:  DRDAR-CLB-PA 
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USER  EVALUATION  OF  REPORT 


Please  take  a  few  minutes  to  answer  the  questions  below;  tear  out 
this  sheet,  fold  as  indicated,  staple  or  tape  closed,  and  place 
in  the  mail.  Your  comments  will  provide  us  with  information  for 
improving  future  reports. 

1 .  BRL  Report  Number _ 

2.  Does  this  report  satisfy  a  need?  (Comment  on  purpose,  related 
project,  or  other  area  of  interest  for  which  report  will  be  used.) 


3.  How,  specifically,  is  the  report  being  used?  (Information 
source,  design  data  or  procedure,  management  procedure,  source  of 
ideas,  etc.) _ _ 


4.  Has  the  information  in  this  report  led  to  any  quantitative 
savings  as  far  as  man-hours/contract  dollars  saved,  operating  costs 
avoided,  efficiencies  achieved,  etc.?  If  so,  please  elaborate. 


5.  General  Comments  (Indicate  what  you  think  should  be  changed  to 
make  this  report  and  future  reports  of  this  type  more  responsive 
to  your  needs,  more  usable,  improve  readability,  etc.) _ 


6.  If  you  would  like  to  be  contacted  by  the  personnel  who  prepared 
this  report  to  raise  specific  questions  or  discuss  the  topic, 
please  fill  in  the  following  information. 

Name: 


Telephone  Number: 
Organization  Address: 


